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Abstract. This paper is an attempt to classify finite-time singularities of PDEs. 

Q^', Most of the problems considered describe free-surface flows, which are easily observed 

.^fH ' experimentally. We consider problems where the singularity occurs at a point, and 



where typical scales of the solution shrink to zero as the singularity is approached. 
Upon a similarity transformation, exact self-similar behaviour is mapped to the fixed 
point of a infinite dimensional dynamical system representing the original dynamics. 
We show that the dynamics close to the fixed point is a useful way classifying the 
structure of the singularity. Specifically, we consider various types of stable and 
unstable fixed points, centre-manifold dynamics, limit cycles, and chaotic dynamics. 
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Figure 1. A drop of Glycerin dripping througli PolyDimetliylSiloxane (PDMS) near 
snap-ofF [T|. The nozzle diameter is 0.48 cm. 



1. Introduction 



A non-linear partial differential equation (PDE), starting from smooth initial data, will 
in general not remain smooth for all times. Consider for example the physical case 
shown in Fig. [H which we will treat in section [3] below. Shown is a snapshot of one 
viscous fluid dripping into another fluid, close to the point where a drop of the inner 
fluid pinches off. This process is driven by surface tension, which tries to minimise the 
surface area between the two fluids. At a particular point xo,to in space and time, 
the local radius h{z, t) of the fluid neck goes to zero. This point is a singularity of 
the underlying equation of motion for the one-dimensional profile h{z, t). In particular, 
typical length scales near the pinch point go to zero (the minimum radius, at the very 
least, but also the solution's axial extend). This absence of a characteristic scale near 
the singularity is the basic motivation to look for self-similar solutions. 

A fascinating aspect of the study of singularities is that they describe a great 
variety of phenomena which appear in the natural sciences and beyond [2]. For 
example, such singular events occur in free-surface flows [3J, turbulence and Euler 
dynamics (singularities of vortex tubes [U [5] and sheets [6]), elasticity [7], Bose- 
Einstein condensates [8] , non-linear wave physics [9] , bacterial growth [TOl [11] , black- hole 
cosmology fi2[ [13] , and financial markets [H] . 

In this paper we consider equations 

ht = F[h], (1) 

where F[h] represents some (nonlinear) differential or integral operator. For simplicity, 
we will for the most part discuss the case of scalar h, but make reference to the many 
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important cases where the variable h is a vector, or ([T]) is a system of equations. Let us 
suppose that ([1]) forms a locahsed singularity at xq, to- If t' = to ~ ^ and x' = x — Xq, we 
are looking for local solutions of ([T]) which have the structure 

h{x,t)=t'"(j){x'/t"'), (2) 

with appropriately chosen values of the exponents a,p. 

Giga and Kohn fibl [T6] proposed to introduce self-similar variables r = — In(t') 
and C, = x' /t'^ to study the asymptotics of blow up. Namely, putting 

h{x,t)=rH{^,r), (3) 

([T]) is turned into the "dynamical system" 

H, = G[H] = aH + t'^-"F[t"'H]. (4) 

If ([2]) is indeed a solution of ([T]), the right hand side of (JH) is independent of r, and 
self-similar solutions of the form ([2]) are fixed points of (jl]). By studying the "long-time" 
(r — »• cxo) behaviour of solutions of (jl]) one can study the behaviour near blow-up. For 
the relation of singular PDE problems to the renormalisation group, developed in the 
context of critical phenomena, see [HI [18], for a more computational perspective, see 

m- 

Solutions to the original PDE ([1]) for given initial data can be viewed as orbits in 
some infinite dimensional phase phase, for instance, L^. If the fixed point is an attractor, 
blow-up will be self-similar for some class of initial initial conditions. However, other 
types of attractors (u;-limit sets in the notation which is customary in the context of 
partial differential equations, see [20] and references therein) are frequently observed, 
furnishing a very fruitful means of classifying singularities. In this paper we discuss the 
following cases: 

(i) Stable fixed points 

Initial conditions in some neighbourhood of the fixed point become attracted to it, 
so the solution converges exponentially to the self-similar solution ([2]). Formally, 
two eigenvalues of the linearisation around the fixed point are always positive, but 
these unstable motions can be absorbed into a redefinition of the origin of space and 
time; this is discussed in section 12.11 A sub-classification into self-similarity of the 
first and of the second kind is due to Barenblatt [2T]. Self-similar solutions are of 
the first kind if ([2]) only solves ([1]) for one set of exponents a, f3; their values are fixed 
by either dimensional analysis or symmetry, and are thus rational. Solutions are of 
the second kind (in the sense of Barenblatt) if solutions ([2]) exist for a continuous 
set of exponents a, (3; the exponents are fixed by a non-linear eigenvalue problem, 
and take irrational values in general. 

(ii) Unstable fixed points 

The linearisation around the fixed point possesses positive eigenvalues, so it is never 
reached, except for non-generic initial data. Often a stable fixed point is associated 
with an infinite sequence of unstable fixed points. 
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(iii) Travelling waves 

Solutions of ([1]) converge to h = t'°'(l){^ + ct), which is a travelhng wave solution of 
(jlj) with propagation velocity c. 

(iv) Centre manifold. 

This is also known, when it leads to singularities that develop at a faster rate than 
the selfsimilar scaling, as type-II self-similarity [22]; it arises if one eigenvalue of the 
linearisation around the fixed point is zero, and there is a non-linear dependence 
instead. This typically leads to corrections involving logarithmic time r. We discuss 
two different cases, involving quadratic and cubic non-linearities. 

(v) Limit cycles 

This is also known as "discrete self-similarity" [121 123]. Corresponding solutions 
have the form h = t'°'ip [^, r] with ip being a periodic function of period T in r. 
Thus at the discrete sequence of times r„ = tq + nT, which approaches the singular 
time for n —>■ oo, the solution looks like a simple self-similar one. 

(vi) Strange attractors 

In principle, more complex behaviour is possible, where the orbits of the dynamical 
system lie on a strange attractor. At the moment, we are not aware of an equation 
exhibiting such behaviour which would correspond to any physical phenomenon. 
However, this may simply be due to the fact that the corresponding singular 
behaviour is more difficult to detect numerically, and to grasp analytically. To 
demonstrate that such behaviour is at least possible, we show that any finite 
dimensional dynamical system may be "embedded" in the singular dynamics. As 
an explicit example, we show that the phase-space trajectory may lie on the Lorenz 
attractor. 

(vii) Exotic objects 

There might be other types of behaviour that have no analogue in finite-dimensional 
dynamical systems. In particular, blow-up may occur at several points (xo,to at 
the same time, in which case the description (jlj) is not so useful. 

This paper's aim is to assemble the body of knowledge on singularities of equations 
of the type ([1]) that is available in both the mathematical and the applied community, 
and to categorise it according to the types given above. In addition to rigorous results 
we pay particular attention to various phenomenological aspects of singularities which 
are often crucial for their appearance in an experiment or a numerical simulation. For 
example, what are the implications of the type of singularity for the approach of the 
PDE solution onto the self-similar form ([2])? In most cases, we rely on known examples 
from the mathematical physics literature. To find an explicit example for limit cycle 
behaviour as well as chaotic dynamics, we propose a new set of model equations, inspired 
by a problem in general relativity. 
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Figure 2. Scanning electron microscopy images illustrating the pinch-ofF of a row of 
rectangular troughs in silicone (top) [24j . The bottom picture shows the same sample 
after 10 minutes of annealing at 1100°C. The troughs have pinched off to form a row 
of almost spherical voids. The dynamics is driven by surface diffusion. 



2. Stable and unstable fixed points 

2.1. Self-similarity of the first kind 

Our example, exhibiting self-similarity of the first kind (in the sense of Barenblatt) 
is that of a solid surface evolving under the action of surface diffusion. Namely, atoms 
migrate along the surface driven by gradients of chemical potential, see Fig|2l The 
resulting equations in the axisymmetric case, where the free surface is described by the 
local neck radius h{x,t), are [25] : 



where 



K 
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h 


h 


[{i + hiy/^ 


1 



Kj'J 



h^ 



(5) 



(6) 



h{l + hiy/^ (1 + /.2)3/2 

is the mean curvature. In ([5]), ([6]), all lengths have been made dimensionless using an 
outer length scale R (such as the initial neck radius), and the time scale R^/D^, where 
D4 is a forth-order diffusion constant. 

At a time t' <C 1 away from breakup, dimensional analysis implies that i = t'^'^ is 
a local length scale. This suggests the similarity form 

h{x,t)=t"^''Hix'/t'^/^), (7) 
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Figure 3. The approach to the self-similar profile for equation ([5]). The dashed line 
is the stable similarity solution H{^) as found from ([8]). The full lines are rescaled 
profiles found from the original dynamics ([5]) at hmin = 10~^, 10~^, and hmin = 10~^, 
respectively. As the singularity is approached, they converge rapidly onto the similarity 
solution d?]). 



implying a = (3 = 1/4. This is the classical situation for self-similarity of the first kind, 
more examples are found in [211 126]. The similarity form of the PDE becomes 



V H 



H 



t'2^ 



-K 



i 



X' 



(8) 



[i + nyi^ 

where the prime denotes differentiation with respect to i^. 

Solutions of ([8]) have been studied extensively in [27j. To ensure matching to a time- 
independent outer solution, the leading order time dependence must drop out from ([7]), 
implying that 

^(0~ciei, e--±oo. (9) 

It turns out that all similarity solutions are symmetric, so only one constant c needs to 
be determined. Exactly as in the closely related problem of surface-tension-driven fluid 
pinch-off, the requirement of a certain growth condition ([9]) is enough to fix a unique 
solution of ([8]) [28j. The value of c comes out as part of the solution. Solutions of ([8]) 
with the growth condition ([9]) form a discretely infinite set [27J, again like the fluids 
problem [29]. The series of similarity solutions is conveniently ordered by descending 
values of the minimum, see table [H 

Next we turn to the dynamical system that describes the dynamics away from the 
fixed point, by putting 

h{x,t)=t''"H{tr), (10) 
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0.701595 


1.03714 


1 


0.636461 


0.29866 


2 


0.456842 


0.18384 


3 


0.404477 


0.13489 


4 


0.355884 


0.10730 


5 


0.326889 


0.08942 



Table 1. A series of similarity solutions of ^ as given in [27j. The higher-order 
solutions become successively thinner and flatter. 



where r = — In(t'). The similarity form of ([5]) becomes 

(11) 



Hr = liH-^H^) + ^ 



H 

'1 + Hiy/^''^ 



which reduces to dH]) if the left hand side is set to zero. To assure matching of (ITTl) to 
the outer solution, we require the boundary condition 

Hr-{H - iH^)/4: -^ for |^| -^ oo. (12) 

Next we linearise (ITTl) around H = H{^), by writing H = H{^) + eP(.^,r), which 
gives 

Pr = C(H)P. (13) 

Since H satisfies (IT2ll . P{C,,t) must do the same. In particular, this means that if 

C(H)P, = u,P,, (14) 

i.e. if z/j is an eigenvalue of the linear operator, the corresponding eigenfunction must 
grow like 

P^i0^e-"''■ (15) 

If the similarity solution H{^) is to be stable, the eigenvalues of C{H) must be 
negative. However, there are always two positive eigenvalues, which are related to the 
invariance of the equation of motion ([5]) under translations in space and time. Namely, 
for any e, the translated similarity solution 

hix,t) = H{^) (16) 

is an equally good self-similar solution of (jS]), and thus of fITT]) . In particular, we can 
expand (fT6l) to lowest order in e, and find that 

H,{^, r) = H{0 + ee^-H'iO + 0{e') (17) 

is a solution of flT^ . 

Thus, since CH = 0, 

ee^r0 ^ dHM, r) ^ ^^^ ^ ee^^£F'. (18) 

OT 
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But this means that u^^ = /3 = 1/4 is an eigenvalue of C with eigenfunction H (^). 
Similarly, considering the transformation t —* t+e, one finds a second positive eigenvalue 
z/f = 1, with eigenfunction $^H . To reiterate, the physical meaning of these eigenvalues is 
that upon perturbing the similarity solution, the singularity time as well as the position 
of the singularity will change. Thus if the coordinate system is not adjusted accordingly, 
it looks as if the solution would flow away from the fixed point. If, on the other hand, 
the solution is represented relative to the perturbed values of xq and to, the dynamics 
will converge onto a stable similarity solution. 

The eigenvalues of the solutions Hi have been found numerically in [27]. The 
result is that the linearisation around the "ground state" solution Hq only has negative 
eigenvalues (apart from the two trivial ones), while all the other solutions have at least 
one other positive eigenvalue. This means that Hq is the only similarity solution that 
can be observed, all other solutions are unstable. Close to the fixed point, the approach 
to Hq will be dominated by the largest negative eigenvalue ui: 

h{x,t)=t"/'[H{0 + ef'^P,{0]. (19) 

For large arguments, the point C,cr where the correction becomes comparable to the 
similarity solution is .^ ~ et"^'^^^~^'^'^, and thus $,cr ~ t'^^^. This means that the region of 
validity of H{C,) expands in similarity variables, and is constant in real space. This rapid 
convergence is reflected by the numerical results reported in Fig. [31 More formally, one 
can say that for any e there is a 5 such that 

\h{x,t)-t"/U{0\<^ (20) 

if |x'| < (5 uniformly as t' -^ 0. 

2.2. Self- similarity of the second kind 

In the example of the previous subsection, the exponents can be determined by 
dimensional analysis, and therefore assume rational values. As Barenblatt |2I] points 
out, there are problems where the scaling behaviour depends on external parameters, 
set for example by the initial conditions. In that case, the scaling exponent can assume 
any value. Often, this value is fixed by some intrinsic property of the equation, resulting 
in an irrational answer. We will call this situation self-similarity of the second kind (in 
the sense of Barenblatt). A particularly simple example of this kind of singularity is 
the pinch-off of a very viscous thread of liquid [301 [3], which we present now. Another 
recent example is the pinch-off of a two-dimensional inviscid sheet |31j . 

For simplicity, we confine ourselves to the case of a slender viscous filament without 
inertia, for which the equation becomes: 

The typical velocity scale 7/77, where 7 is the surface tension and rj is the viscosity, 
has been absorbed into the time variable. The particularly simple form of (!2T|) has 
been achieved by writing the thread radius in Lagrangian variables, i.e. as function of a 
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Figure 4. A drop of viscous fluid falling from a pipette 1 mm in diameter [33] . Note 
the long neck. 



particle label s. This means the particle is at position z{s, t) at time t, and Zt{s, t) is the 
velocity at time t. The time-dependent constant of integration C{t) must be determined 
from the constraint that u{s,t) = Zg = l/h'^{s,t). 

Note that the self-similar form ([2]) is a solution of ( 1211) for a = 1, and any value of 
(3. Thus contrary to the example described in the previous section, the exponents are 
not completely determined from dimensional analysis or from balancing powers of t' in 
the equation of motion. This is a typical situation in which self-similarity of the second 
kind is observed [32]. Instead, the unknown exponent is determined from a non- linear 
eigenvalue equation, and takes an irrational value. 

Since a = 1 we introduce 

u = t'-^f{^), with i = s/t'^ (22) 



and 



C{t) = Kt' . (23) 



Hence 

where K is an arbitrary constant. Imposing symmetry and regularity of /, we introduce 
an expansion of f{S,) of the form 

/(O = i?o"' + Ce'" + O(e'"), n = l,2,... 

into ( 12^ to obtain the condition 

and define (3 = n(3. Equation f[24l) can easily be integrated in terms of In^ and y = ^/f: 

'^^ ^ lne + C = ^lne" + C, 



((l + 6i?o)2/=^-y2-6i?oy) 6i?o/3 6i?o/? 

with C an arbitrary constant. Computing the integral above we obtain 

y-P [[2^-l)y+lf-^^ {l-yY = CC, (26) 
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where we can fix, without loss of generality, C = 1. In terms of the profile h{s,t) this 
corresponds to fixing the scale of the spatial variable s. Solutions are undetermined up 
to such a scale factor, as is clear from the invariance of ( 12T|) under a change in spatial 
scale. As a result, the axial scale is fixed by the initial conditions. 
The value of the velocity C„ at infinity is therefore given by 

From the condition that the velocity at infinity must vanish we thus obtain Cn = 0, 
which will be the equation that determines the exponent f5. Taking the derivative of 
( l26l) we obtain 

"«""! = 55 (^"'«'''-')^ + i)'"<i-^'' 



and hence 



Vli^) 



-"(^-ii^i^^-'KGff^^^'-^' 



^+/3 //„^ 



1 2n-2l3 + l 



y-—{{2P-l)y + l) 5 " {l-yr-^^^]dy. (27) 

The function Kn{(3) may be written explicitly as 

TiA-/3)T{^)2nP-l 



Kn{i3) = P- 



r(4-/?+^) 2n/?-2 



2r7 -I- 1 1 

—^-P,4-p-4-P+—-l-2nP 

^T{3-P+^) V 2n '^' '^' '^ 2ri' "^ J ^ ^ 

with roots given by table [21 If one converts the Lagrangian variables back to the original 
spatial variables, one obtains 

h{x,t)=t'<Pst{x'/t'^-'). (29) 

Thus for t' ^ the typical radial scale t' of the generic n = 1 solution rapidly becomes 
smaller than the axial scale t"^-^'^^ (cf. table [2]). This explains the long necks seen in 
Fig.H 

For generic initial data Uq{s) = hQ^{s) = Bq + Bis"^ + i?2S'' + ... + Bjs"^^ + ... one 
expects that Bi ^ 0, so that the self-similar solution with n = \ will develop. Only 
if Sj = for i = 1,2, ...,n — 1 and i?„ 7^ the n-th self-similar solution will be the 
asymptotic description of the solution. For this reason, only the n = \ solution is 
stable, since a generic perturbation of the initial data with Bi = for i = l,2, ...,n — 1 
will, in general, make Bi ^ 0. 
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n (3 



1 


2.1748 


2 


2.0454 


3 


2.0194 


4 


2.0105 


5 


2.0065 


10 


2.0014 



Table 2. A list of exponents, found from Kn{f3) = 0, with Kn given by (I28p . 
The number 2n gives the smallest non- vanishing power in a series expansion of the 
corresponding similarity solution around the origin. Only the solution with n = 1 is 
stable. The minimum radius is found from ((25|) . 



3. Travelling wave 

The pinching of a liquid thread in the presence of an external fluid is described by the 
Stokes equation [31] . For simplicity, we consider the case that the viscosity t] of the fluid 
in the drop and that of the external fluid are the same. An experimental photograph 
of this situation is shown in Fig. [H To further simplify the problem, we make the 
assumption (the full problem is completely analogous) that the fluid thread is slender. 
Then the equations given in [1] simplify to 

ht = -v,h/2 - vK, (30) 

where 

''~4A_ mz') + {z - z'yy/^ ^ ' 

Here we have written the velocity in units of the capillary speed Vjj = 'j/v- The limits 
of integration Z- and z+ are for example the positions of the plates which hold a liquid 
bridge [35] . 

Dimensionally, one would once more expect a local solution of the form 

hiz,t)=t'H,^Jj\ (32) 

and HoutiO ^^^ ^'^ be a linear function at infinity to match to a time- independent outer 
solution. In similarity variables, (I5T1) has the form 

v;,t(0 = 7 / °''' ^dt (33) 

We have chosen Zb as a real-space variable close to the pinch-point, such that the 
similarity description is valid in [— Zb,^^]. But if i^out is linear, the integral in (l33ll 
diverges, which means that a simple "fixed point" solution (l32l) is impossible. 

However, the integral can be made convergent by introducing a shift in the similarity 
variable ^: 

/i = t'ifout(e-6r), (34) 
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with T = — In(t') as usual. This means in similarity variables the solution is a travelling 
wave. With this modification, the mass balance fl30l) becomes 

- //out + KA^ + Vont + hr) = HouXJ2. (35) 

Now we choose b such that the logarithmic singularity cancels, namely we demand that 

' ^ouA^) d^'-b\n{t') (36) 



(37) 



4 J-z,/t' ^/Hl, + (e - e)2 
finite for t' —>■ 0. This is achieved by putting 



Thus defining 



H, H 



'HI + 1 JHl + 1 



Vfin(0= lim - [ ^°"' ^d^' + b\nA (38) 

the similarity equation 

- ifo,t + i^;ut(e + ^fin - eo) = ^outKn/2 (39) 

is finite, and ^o is an arbitrary constant. It remains as an arbitrary axial shift in the 
similarity solution. 

The numerical solution of the integro-differential equation (l39l) gives 

hmin = CLontVrit', whcrC Oout = 0.0335. (40) 

The slope of the solution away from the pinch-point are given by 

H+ = 4.81 and H^ = -0.105, (41) 

which means the solution is very asymmetric, as confirmed directly from Fig. [H 

4. Centre manifold 

In section |2] we described the generic situation that the behaviour of a similarity solution 
is determined by the linearisation around it. In the case of a stable fixed point, 
convergence is fast, and the observed behaviour is essentially that of the fixed point. 
In this section, we describe two different cases where the largest eigenvalue vanishes, so 
higher-order non-linear terms have to be taken into account. The approach to the fixed 
point is now much slower, and much more of the observed behaviour is determined by 
the approach to the fixed point. Depending on the type of non-linearity, there is great 
freedom of possible behaviours, of which we discuss two. 
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Figure 5. Nine images (of width 3.5 mm) showing how a '^He crystal "flows" down 
from the upper part of a cryogenic cell into its lower part [37) . The recording takes a 
few minutes, the temperature is 0.32 K. 11 mK. The crystal first "drips" down, so that 
a crystalline "drop" forms at the bottom (a to c); then a second drop appears (d) and 
comes into contact with the first one (e); coalescence is observed (f) and subsequently 
breakup occurs (h). 



4.1. Quadratic non-linearity 

Axisymmetric motion by mean curvature in three spatial dimensions is described by the 
equation 

i^xx t 

where h{x,t) is the radius of the moving free surface. A very good physical realization 



ht 



(42) 



of 



is the melting and freezing of a '^He crystal, driven by surface tension 



see 

Fig. [5l As before, the time scale t has been chosen such that the diffusion constant, 
which sets the rate of motion, is normalised to one. A possible boundary condition for 
the problem is that h{0, t) = h{L, t) = R, where R is some prescribed radius. For certain 
initial conditions h{x,0) = ho{x) the interface will become singular at some time to, at 
which h{xo,to) = and the curvature blows up. The moment of blow-up is shown in 
panel h of Fig. [5l for example. 

Inserting the self-similar solution ([2]) into (l42l) . one finds a balance for a = (3 = 1/2. 
The corresponding similarity equation is 

One solution of (H3ll is the constant solution (/)(,^) = v2- Another potential solution 
is one that grows linearly at infinity, to ensure matching onto a time-independent 
outer solution. However, it can be shown that no solution to fH3l) . which also grows 
linearly at infinity, exists [3S1 EH] • Our analysis below follows the rigorous work in [22] , 
demonstrating type-II self-similarity. In addition, we now show how the description of 
the dynamical system can be carried out to arbitrary order. 
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The solution relevant to our analysis is the constant solution, but which of course 
does not match onto a time-independent outer solution. We thus write the solution as 

h{x,t)=t''/'[V2 + g{^,T)], (44) 

with T = — In(t') as usual. The equation for g is then 

which we solve by expanding into eigenfunctions of the linear part of the operator 

Cg = g-^g'/2 + g". (46) 

It is easily confirmed that 

CH2M/2) = {l-z)H2^{^/2), (47) 

where Hn is the n-th Hermite polynomial [40j: 

Hn = {-iTe- ^e- . (48) 

Thus the linear part of (l45l) becomes 

^ = (l-^)a., (49) 

which means that all eigenvalues are negative except for the first, which vanishes. To 
investigate the approach of the cylindrical solution, one must therefore include nonlinear 
terms in the equation for Oi. 
If we write 

(7(e,r) = J]a,(r)i72.(e/2), (50) 

the equation for ai becomes 

^ = -23/^a? + 0(aia,), (51) 

whose solution is 

ai = l/(23/V). (52) 

Thus instead of the expected exponential convergence onto the fixed point, the 
approach is only algebraic. Since all other eigenvalues are negative, the r-dependence 
of the tti is effectively determined by ai. Namely, as we will see below, aj = 0{t~^), so 
corrections to flSTl) are of higher order. 

If one hnearises around fl52l) . putting ai = a[' + ei, one finds 

— — = Ci + other terms. (53) 

dr T 

This means that the coefficient A oi ei = AJt^ remains undetermined, and a simple 

expansion of Oj in powers of r^^ yields an indeterminate system. Instead, at quadratic 

order, a term of the form e\ = AXw tJt^ is needed. Fortunately, this is the only place in 
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the system of nonlinear equations for Oj where such an indeterminacy occurs. Thus all 
logarithmic dependencies can be traced, leading to the general ansatz 



k=i+l 1=0 



5+EE^*- (54) 



a. 



where bi and 5;^^ are coefficients to be determined. The index n is the order of the 
truncation. 

Indeed, the coefficients can be found recursively by considering terms of successively 
higher order in r~^ in the first equation: 

-^ = -l^l^al - 24:V2aia2 + 22a? - 
dr 

T11\Pia\ - \^\\Pia\ + 192a^a2 (55) 

^ = -aa - \p2l^a\ + 6a? - 8V2aia2 (56) 

(IT 

The next two orders will involve the next coefficient as. From fl55l) and fl56l) . one 
first finds ^121 and ^2, by considering Oir^'^^ and 0{t^'^\ respectively. Then, at order 
Qj^^-(n+i)-j jj-^ ^]^g gj,g^ equation, where n = 3, one finds all remaining coefficients bi^i 
in the expansion (!54l) up to fc = n. At each order in r~^, there is of course a series 
expansion in In r which determines all the coefficients. 

We constructed a MAPLE program to compute all the coefficients up to arbitrarily 
high order (10th, say). Up to third order in r~^ the result is: 

, ^f2 17 In (r) v^ 73 ^2 

ai = l/4 — + -^ - — ^ 

T 16 r^ 16 T-^ 

867 1n(r)v^ 289 (In (r))' v^ 



«^ = -V32^ + ^^-^^i^, (5^ 



128 r3 128 r3 

v^ , 5 V2 17 In (r) v^ 
16 7^~ 64 ^ 
and thus h{x,t) becomes 

h{x, t) = t"/^ [V2 + ai(r) (-2 + e) + a2(r) (l2 - 12^' + e')] , (59) 

from which one finds the minimum. To second order, the result is 



'"rmn 



(2t 



'U/2 



1 3 + 171nr 



(60) 



2r 8r2 

Two remarks are in order. First, the presence of logarithms implies that there is 
some dependence on initial conditions built into the description. The reason is that the 
argument inside the logarithm needs to be non-dimensionalised using some "external" 
time scale. More formally, any change in time scale t = t/tg leads to an identical equation 
if also lengths are rescaled according to h = h/ ^JIq. This leaves the pref actor in fl60l) 
invariant, but adds an arbitrary constant tq to r. This is illustrated by comparing to a 
numerical simulation of the mean curvature equation fH21) close to the point of breakup, 
see Fig. [61 Namely, we subtract the analytical result (!60|) from the numerical solution 
hmin/{'^V^) and multiply by r^. As seen in FigJHl the remainder is varying slowly over 12 
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Figure 6. A plot of 



To)/8) (full line) with To = 4.56. 



h^^n/V^t' - 1 + 1/(2t) t2 (dashed line) and ro/2- (3+17 ln(T- 



decades in t' . If the constant tq is adjusted, this small variation is seen to be consistent 
with the logarithmic dependence predicted by ( 160|) . 

The second important point is that convergence in space is no longer uniform as 
implied by (120|) for the case of self-similarity of the first kind. Namely, to leading 
order the pinching solution is a cylinder. For this to be a good approximation, one 
has to require that the correction is small: ■C^/r ^ 1. Thus corrections become 
important beyond ^cr ~ t, which, in view of the logarithmic growth of r, implies 
convergence in a constant region in similarity variables only. As shown in [36], the 
slow convergence toward the self-similar behaviour has important consequences for a 
comparison to experimental data. 



4-2. Cubic non-linearity 

The next example is that of bubble breakup [H], for which a very different form of 
nonlinearity is observed. As shown in [H], the equation for a slender cavity or bubble 
is 



a(e,t)rfe 



a" 



(61) 



_L^[z-iy + a{z,t) 2a' 

where a{z,t) = h'^{z,t). The integral runs over the fluid domain. If for the moment 
one disregards boundary conditions looks for solutions to (16T]) of cylindrical form, 
a{z,t) = ao{t), one can do the integral to find 



Ooln 



4L2 



a2 



2ao' 



(62) 
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Figure 7. The pinch-ofF of an air bubble in water [44j. An initially smooth shape 
develops a localised pinch-point. 




Figure 8. A comparison of the exponent a between full numerical simulations of 
bubble pinch-off (solid line) and the leading order asymptotic theory (p5)) (dashed 
line). 



It is easy to show that an an asymptotic solution of fl62|) is given by 

At . . 

"° °^ MA^' (^^) 

corresponding to a power law with a small logarithmic correction. Indeed, initial theories 

of bubble pinch-off [l2l HS] treated the case of an approximately cylindrical cavity, which 

leads to the radial exponent a = 1/2, with logarithmic corrections. 

However both experiment [H] and simulation [51] show that the cylindrical solution 
is unstable; rather, the pinch region is rather localised, see Fig. [71 Therefore, it is not 
enough to treat the width of the cavity as a constant L; the width A is itself a time- 
dependent quantity. In [41j we show that to leading order the time evolution of the 
integral equation ( 1611) can be reduced to a set of ordinary differential equations for the 
minimum oq of a{z,t), as well as its curvature Oq. 

Namely, the integral in ( 1611) is dominated by a local contribution from the pinch 
region. To estimate this contribution, it is sufficient to expand the profile around the 
minimum at z = 0: a{z,t) = qq + aQ/2z^ -|- 0{z^). As in previous theories, the integral 
depends logarithmically on a, but the axial length scale is provided by the inverse 
curvature A = (2ao/ao)^^^- Thus evaluating (I6T1) at the minimum, one obtains [H] to 
leading order 

aoln(4AVao) = ag/(2ao), (64) 
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which is a coupled equation for oq and A. Thus, a second equation is needed to close the 
system, which is obtained by evaluating the the second derivative of fl6Tl) at the pinch 
point: 

«oln (4t) _ 2^ = ^ - M<. (65) 

\e^a'^J ao qq 2ag 

The two coupled equations fl64l) . fl65|) are most easily recast in terms of the time- 
dependent exponents 

2a = —drao/ao, 26 = —dra'^/aQ, (66) 

where r = — Int' and (3 = a — 6, are a generalisation of the usual exponents a and (3. 
The exponent S characterises the time dependence of the aspect ratio A. Returning to 
the collapse (|62|) predicted for a constant solution, one finds that a = 1/2 and 6 = 0. 
In the spirit of the the previous subsection, this is the fixed point corresponding to 
the cylindrical solution. Now we expand the values of a and S around their expected 
asymptotic values 1/2 and 0: 

a = 1/2 + u{t), S = v{t). (67) 

To leading order in 6, the resulting equations are 

drU = -Svu"^, drV = Sv^ , (68) 

which describe perturbations around the leading-order similarity solution. These 
equations are analogous to (ISTI) . but they have a degeneracy of third order, rather 
than second order. Equations (l68ll are easily solved to yield, in an expansion for small 

a = 1/2 + ^ + 0(r), 6 = ^ + 0{r-y'). (69) 

Thus the exponents converge toward their asymptotic values a = (3 = 1/2 only very 
slowly, as illustrated in Fig. [HI This explains why typical experimental values are found 
in the range a ~ 0.54 — 0.58 [^, and why there is a weak dependence on initial 
conditions [15]. 

The cubic equation (1671) applies to the exponents a, 6, rather than the solution itself, 
as in the previous subsection, where we dealt with mean curvature flow. In fact, it is 
easy to re-analyse the solution to the mean curvature problem, and to formulate it in 
terms of time-dependent exponents. Let us define Oq and Oq for the mean curvature 
problem through h'^{z,t) = ao + ag/2z^ + O(z^). From ( l44l) it follows that ao = 2t' and 
a'o = 2(^2 + g)g". This gives a^ = A^/2ai, and thus 

dj = -26\ (70) 

instead of ( 168|) . This means mean curvature flow, formulated in terms of exponents, once 
more gives a quadratic non-linearity, rather than the cubic term (168|) found for bubble 
collapse. From (jTOjl one finds to leading order 6mc = 1/(2t) and a^c = 1/2 + 0{1/t'^) 
for the mean curvature flow. 
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5. Limit cycles 

An example for this kind of blow-up was introduced into the literature in [1^ in the 
context of cosmology. There is considerable numerical evidence |16] that discrete self- 
similarity occurs at the mass threshold for the formation of a black hole. The same 
type of self-similarity has also been proposed for singularities of the Euler equation [17] 
and for a variety of other phenomena [H] . A reformulation of the original cosmological 
problem leads to the following system: 

A ^ ^^, (71) 

ia-% ^ ^-'^+^^ + "-^V''\ (72) 



X 



ia-\ 



if + x)U' - if - x)V 



2 



+ 1 



X 



A' - 1, (73) 



xif +x) 
y^^ mi-a^W.VH.V. ^ (75) 

xif - x) 

In [13], the self-similar description corresponding to the system (I7T1) -( 1751) was solved 
using formal asymptotics and numerical shooting procedures. This leads to the solutions 
observed in p!2]. Below we propose a very similar system, which we solve analytically, 
and which shows the same type of limit cycle behaviour: 



ut^u^ = 2fv, (76) 

vt + v^ = -2fu, (77) 

ft = f. (78) 

The simplified system ( 176l) -( l78l) can be solved introducing characteristics: 

1] = X + t, u = X — t, (79) 

which leads to 

Uri = fv, V^ = -fu. (80) 

The system (1801) has to be solved at constant v with initial conditions 

m(z/, v) = Uinitil^), Viu, V) = Vinitil^), (81) 

where Umitix) = uix, 0) and Vmuix) = f (x, 0) are the profiles at t = 0. 
The solution of (1751) is 



^<-"^l77w^' <''' 
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where /j„it(x) is the initial profile, which can be written in terms of u,-)] as 

1 



fiv,^) 



l//.™t((r7 + z/)/2)-(r7-z/)/2- 
The function / can be eliminated from flHOl) using the transformation 



J V 

which results in the system for a harmonic oscillator: 

Uq = V, V(^ = —U. 

The solution is 

u = C{u) sin(C + (pQ^)), V = C{iy) cos(C + M^^)), 

with C(z/) = ^/u^niM +ViniM and (f)oiu) = arcsm^Uinuii^) / C (u)) . 

According to (!82|) . a singularity of the PDE system (!76l) - (178l) first occurs at the 
maximum /o of finit{x), which we assume to occur at x = without loss of generality. 
Thus locally we can write f{x, 0) ~ /o — ax^, and for small t' = to — t we have 

t'-i 



(83) 

(84) 

(85) 
(86) 



/ 



e 



X 



(87) 



1 + atle' t''/'' 

for any finite ^. Using this explicit form of /, (!84|) can be integrated to find C- At 
constant ^ we have 

iy = ~to + it'^'^ + t', r/ = to + it'^'^ - t', (88) 

so taking the limit t' — > 0, the leading order result for f l84l) is 



c 



In 



t'{l + atl^'^){l + at, 



T+m- 



(89) 



Thus as the singularity is reached, t' — > 0, the variable ( goes to infinity. This means 
the singularity corresponds to the long-time limit of the dynamical system fl85l) . which 
will perform a harmonic motion according to ( l86l) . Namely, for t' ^ the result is 



u 



'^init\ 



to) + vf„J-to) sin < arcsin 



Uinity—to) 



V'^lniti-to) + Vfniti-to) 



In 



t^l + at^e )(1 + Qt 
to 



Thus the singular solution is of the general form 



(90) 



(91) 



where ip is periodic in r. This is a particularly simple version of discretely self-similar 
behaviour. Note that the character of the solution is different from the travelling waves 
described in section [31 
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6. Strange attractors and exotic behaviour 

In connection to limit cycles and in the context of singularities in relativity, a few 
interesting situations have been found numerically quite recently. One of them is the 
existence of Hopf bifurcations where a self-similar solution (a stable fixed point) is 
transformed into a discrete self-similar solution (limit cycle) as a certain parameter 
varies (see [19]). Other kinds of bifurcations, for example of the Shilnikov type, are 
found as well [50]. Now we demonstrate that chaotic behaviour is also possible. 

In section 14.11 we treated a system of an infinite number of ordinary differential 
equations for the coefficients of the expansion of an arbitrary perturbation to an explicit 
solution. Such high-dimensional systems in principle allow for a rich variety of dynamical 
behaviours, including those found in classical finite dimensional dynamical systems, such 
as chaos. Consider for instance an equation for the perturbation g (the analogue of (HSj)) 
of the form 

gr = Lg + F{g,g), (92) 

where Lg is a linear operator. Assuming an appropriate non-linear structure for the 
function F, an arbitrary nonlinear (chaotic) dynamics can be added. 

To give an explicit example of a system of PDEs exhibiting chaotic dynamics, 
consider the structure of the example given in the previous section [51 It can be 
generalised to produce any low-dimensional dynamics near the singularity. Namely, 
let us generaUse the system (!76l) - (!78l) to 

w«+«« = 2/F,({««}), 2 = l,...,n, (93) 

/* = f. (94) 

Using the transformation ( 1791) . the first n equations are turned into: 

uf/f = F,, ^ = l,...,n, (95) 

which is an ODE system for constant v. The system (l95l) has to be solved with initial 
conditions 

w«(z.,z/)=4),(z.), z = l,...,n. (96) 

As before, the function / can be eliminated using (18^ . resulting in the general 
non-linear dynamical system 

uf = Fi, 2 = l,...,n. (97) 

Now if one chooses n = 3 and 

F, = a(n(2)-n«), F2 = pu^^^ -u^'^ -u^'^ u'^'l F3 = u^'^'^'^ - (3u^^\ (98) 



(I97j) is the Lorenz system [51] . 

As before, for t' —>■ 0, the variable ( goes to infinity, and near the singularity one 
is exploring the long-time behaviour of the dynamical system (197]) . In the case of fl98|) . 
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and for sufficiently large p, the resulting dynamics will be chaotic. Specifically, taking 
cr = 10, p = 28, and (3 = 8/3, as done by Lorenz [52], the maximal Lyapunov exponent is 
0.906. Now if U^^\u, (), f/(2)(z/, (), U^^\u, () is a solution of ([97]) with initial conditions 
(l96|l . the final form of the similarity solution is 

M«(a;^t) = t/«(-to + e^"/',r + 0(O), 2 = 1,2,3. (99) 

In the limit t' ^ one cannot replace the first argument by —to + ^t'^^'^ ~ —to- The 
reason is that in the limit of large r two trajectories diverge like exp(AmaxT) = t'^'^™""', 
where the largest Lyapunov exponent is larger than 1/2. 

7. Outlook 

The singularities described in this paper are point-like in the sense that they occur at 
a single point xq at a given time to. There are two important situations we have not 
discussed and for which the dynamical systems point of view and analytical approach 
does not apply. First, the case of singularities not developing at single points, but on 
sets of finite measure. This is the case in a few simple examples of reaction-diffusion 
equations of the family 

Ut-Au = uP-b \Vu\'^ for x eQ (100) 

where depending on values of p > 1 and g > 1 singularities in the form of blowing-up u 
may be regional {u blows up in subsets of Q of finite measure) or even global (the solution 
blows-up in the whole domain). See for instance [53] and references therein. Another 
interesting possibility is that the stable fixed point which is approached depends on the 
initial conditions. Thus exponents could vary either discretely or continuously with the 
choice of initial conditions. The infinite sequence of similarity solutions found for ( 12T11 
is not an example for such behaviour, since only one solution is stable. All other fixed 
points will not be visible in practise. 

Singularities may even happen in sets of fractional Hausdorff dimension, i.e., 
fractals. This is the case of the inviscid one-dimensional system for jet breakup (cf. [51] ) 
and might be case of Navier-Stokes system in three dimensions, where the dimension of 
the singular set at the time of first blow-up is at most 1 (cf. [55]). This connects with the 
second issue we did not address here. It is the nature of the singular sets both in space 
and time. In many instances, existence of global in time (for all < t < cxo) solutions to 
nonlinear problems can be established in a weak sense, that is, allowing certain kind of 
singularities to develop both in space and time. In the case of 3-D Navier-Stokes system, 
the impossibility of singularities "moving" in time, that is of curves x = (f{t) in the 
singular set is well-known [55] . Hence, provided a certain singularity does not persist in 
time, the question is how to continue the solutions after a singularity has developed. 
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